#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math
from scipy.stats import skew


# In[2]:


df = pd.read_csv("D:/research/recordbreak/record45/chenab_record45_dec.csv")


# In[3]:


# Assuming 'year' is the name of the column containing the years
df['year'] = pd.to_datetime(df['year'], format='%Y').dt.year
df.set_index('year', inplace=True)


# In[4]:


print(df)


# In[5]:


dr_theo = df[df.index <= 2099][['dr_theo']]
print(dr_theo)


# In[6]:


fr_theo = df[df.index <= 2099][['fr_theo']]
print(fr_theo)


# # Decreasing

# In[7]:


dr_obs_dec = df[df.index <= 2022][['dr_obs']]
print(dr_obs_dec)


# In[8]:


fr_obs_dec = df[df.index <= 2022][['fr_obs']]
print(fr_obs_dec)


# In[9]:


dr_boot_dec = [f"dr_dec{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
dr_boot_dec = df.loc[df.index <= 2099, dr_boot_dec]
print(dr_boot_dec)


# In[10]:


fr_boot_dec = [f"fr_dec{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
fr_boot_dec = df.loc[df.index <= 2099, fr_boot_dec]
print(fr_boot_dec)


# In[11]:


percentile_5_dr = np.percentile(dr_boot_dec, 5, axis=1)
percentile_95_dr = np.percentile(dr_boot_dec, 95, axis=1)
percentile_50_dr = np.percentile(dr_boot_dec, 50, axis=1)


# In[13]:


percentile_5_fr = np.percentile(fr_boot_dec, 5, axis=1)
percentile_50_fr= np.percentile(fr_boot_dec, 50, axis=1)
percentile_95_fr = np.percentile(fr_boot_dec, 95, axis=1)


# In[14]:


# Calculate the ratio between corresponding percentiles of two datasets
ratio_5 = np.divide(percentile_5_fr, percentile_95_dr)
ratio_50 = np.divide(percentile_50_fr, percentile_50_dr)
ratio_95 = np.divide(percentile_95_fr, percentile_5_dr)


# In[15]:


ratio_obs_dec = df[df.index <= 2022][['ratio_obs_dec']]
print(ratio_obs_dec)


# # RCP 8.5

# In[19]:


df1 = pd.read_csv("D:/research/recordbreak/record85/chenab_record85_dec.csv")


# In[20]:


# Assuming 'year' is the name of the column containing the years
df1['year'] = pd.to_datetime(df1['year'], format='%Y').dt.year
df1.set_index('year', inplace=True)
print(df1)


# In[21]:


dr85_boot_dec = [f"dr_dec{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
dr85_boot_dec = df1.loc[df1.index <= 2099, dr85_boot_dec]
print(dr85_boot_dec)


# In[29]:


fr85_boot_dec = [f"fr_dec{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
fr85_boot_dec = df1.loc[df1.index <= 2099, fr85_boot_dec]
print(fr85_boot_dec)


# In[30]:


percentile_5_dr85 = np.percentile(dr85_boot_dec, 5, axis=1)
percentile_95_dr85 = np.percentile(dr85_boot_dec, 95, axis=1)
percentile_50_dr85 = np.percentile(dr85_boot_dec, 50, axis=1)


# In[31]:


percentile_5_fr85 = np.percentile(fr85_boot_dec, 5, axis=1)
percentile_95_fr85 = np.percentile(fr85_boot_dec, 95, axis=1)
percentile_50_fr85 = np.percentile(fr85_boot_dec, 50, axis=1)


# In[32]:


ratio85_5 = np.divide(percentile_5_fr85, percentile_95_dr85)
ratio85_50 = np.divide(percentile_50_fr85, percentile_50_dr85)
ratio85_95 = np.divide(percentile_95_fr85, percentile_5_dr85)


# # New

# In[33]:


plt.figure(figsize=(16, 9))
# Plotting the drought record-breaking events decay rate
plt.plot(dr_obs_dec.index, dr_obs_dec.values, label='Observed', color='black', alpha=0.8, linewidth=4)

# plt.plot(dr85_m_dec.index, dr85_m_dec, label='DNN', color='purple', linewidth=3) #, linestyle='dotted'
#for col in dr_boot_dec.columns:
    #plt.plot(dr_boot_dec.index, dr_boot_dec[col], color='teal', alpha=0.03)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_dec.index, percentile_50_dr, color='teal', linewidth=3, alpha=0.9, label='RCP 4.5')
#plt.plot(dr_boot_dec.index, percentile_5_dr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5)
#plt.plot(dr_boot_dec.index, percentile_95_dr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')

plt.fill_between(dr_boot_dec.index, percentile_5_dr, percentile_95_dr, color='teal', alpha=0.2)

#for col in dr85_boot_dec.columns:
    #plt.plot(dr85_boot_dec.index, dr85_boot_dec[col], color='red', alpha=0.03)

# Plot 5th and 95th percentiles
plt.plot(dr85_boot_dec.index, percentile_50_dr85, color='red', linewidth=3, alpha=0.8, label='RCP8.5')
#plt.plot(dr85_boot_dec.index, percentile_5_dr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5)
#plt.plot(dr85_boot_dec.index, percentile_95_dr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')

plt.fill_between(dr_boot_dec.index, percentile_5_dr85, percentile_95_dr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theoretical (1/n)', color='blue', linewidth=3)
# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=22,fontweight='bold') #fontweight='bold'
plt.ylabel('Probabaility (Log)', fontsize=22,fontweight='bold')
#plt.title('Chenab drying Annual-Min Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=20,fontweight='bold')
plt.yticks(fontsize=18,fontweight='bold')
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=8, width=2)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
plt.legend( prop={'size': 14, 'weight': 'bold'})

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=18)
for label in ax.yaxis.get_ticklabels():
    label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(2)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=8, width=2)  # Set the desired length, e.g., 10

# Display the plot
plt.show()


# In[44]:


#plt.title('Chenab drying Annual-Min Record-breaking', fontsize=20, fontweight='bold')
plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(dr_obs_dec.index, dr_obs_dec.values, label='Observed', color='black', alpha=0.9, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_dec.index, percentile_50_dr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(dr_boot_dec.index, percentile_5_dr, percentile_95_dr, color='teal', alpha=0.2)

plt.plot(dr85_boot_dec.index, percentile_50_dr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(dr_boot_dec.index, percentile_5_dr85, percentile_95_dr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theoretical (1/n)', color='blue', linewidth=3, alpha=0.8)
# Adding labels and title
#plt.xlabel('Water Year(Apr-Mar)', fontsize=22) #fontweight='bold'
#plt.ylabel('Probabaility (Log)', fontsize=22)


# Increase axis ticks size
plt.xticks(fontsize=30) #,fontweight='bold'
plt.yticks(fontsize=30)
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
#plt.legend( prop={'size': 30}) #, 'weight': 'bold'

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
#ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=30)
#for label in ax.yaxis.get_ticklabels():
    #label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=16, width=4)  # Set the desired length, e.g., 10
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[43]:


plt.figure(figsize=(16, 9))
# Plotting the flood record-breaking events decay rate
plt.plot(fr_obs_dec.index, fr_obs_dec.values, label='Observed', color='black', alpha=0.9, linewidth=4)

#for col in fr_boot_dec.columns:
     #plt.plot(fr_boot_dec.index, fr_boot_dec[col], color='teal', alpha=0.05)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_dec.index, percentile_50_fr, color='teal',  alpha=0.9,linewidth=3, label='RCP 4.5')
#plt.plot(fr_boot_dec.index, percentile_5_fr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
#plt.plot(fr_boot_dec.index, percentile_95_fr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5)

plt.fill_between(fr_boot_dec.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

# 100 predictions from bootstrapping
#for col in fr85_boot_dec.columns:
    #plt.plot(fr85_boot_dec.index, fr85_boot_dec[col], color='red', alpha=0.05)

# Plot 5th and 95th percentiles
plt.plot(fr85_boot_dec.index, percentile_50_fr85, color='red',  alpha=0.8,linewidth=3, label='RCP8.5')
#plt.plot(fr85_boot_dec.index, percentile_5_fr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5,label='5-95 CI')
#plt.plot(fr85_boot_dec.index, percentile_95_fr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5)
plt.fill_between(fr_boot_dec.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theoretical (1/n)', color='blue', linewidth=3)
# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=22,fontweight='bold') #fontweight='bold'
plt.ylabel('Probabaility (Log)', fontsize=22,fontweight='bold')
#plt.title('Chenab drying Annual-Max Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=20,fontweight='bold')
plt.yticks(fontsize=18,fontweight='bold')
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=8, width=2)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
plt.legend( prop={'size': 14, 'weight': 'bold'})

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=18)
for label in ax.yaxis.get_ticklabels():
    label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(2)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=8, width=2)  # Set the desired length, e.g., 10

# Display the plot
plt.show()


# In[42]:


#plt.title('Chenab drying Annual-Max Record-breaking', fontsize=20, fontweight='bold')
plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(fr_obs_dec.index, fr_obs_dec.values, label='Observed', color='black', alpha=0.9, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_dec.index, percentile_50_fr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(fr_boot_dec.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

plt.plot(fr85_boot_dec.index, percentile_50_fr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(fr_boot_dec.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)

#Theoratical
plt.plot(fr_theo, label='Theoretical (1/n)', color='blue', linewidth=3, alpha=0.8)
# Adding labels and title
#plt.xlabel('Water Year(Apr-Mar)', fontsize=22) #fontweight='bold'
#plt.ylabel('Probabaility (Log)', fontsize=22)


# Increase axis ticks size
plt.xticks(fontsize=30) #,fontweight='bold'
plt.yticks(fontsize=30)
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
#plt.legend( prop={'size': 30}) #, 'weight': 'bold'

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
#ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=30)
#for label in ax.yaxis.get_ticklabels():
    #label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=16, width=4)  # Set the desired length, e.g., 10
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[46]:


plt.figure(figsize=(20,10))

# Calculate the inverse of the values
#inv_ratio_obs_dec = 1 / ratio_obs_dec.values
inv_ratio_50 = 1 / ratio_50
inv_ratio_5 = 1 / ratio_5
inv_ratio_95 = 1 / ratio_95

inv_ratio85_50 = 1 / ratio85_50
inv_ratio85_5 = 1 / ratio85_5
inv_ratio85_95 = 1 / ratio85_95

# Plotting the inverse of the flood record-breaking events decay rate
plt.plot(ratio_obs_dec.index, ratio_obs_dec, label='Observed', color='black', alpha=0.9, linewidth=5)
plt.plot(fr_boot_dec.index, inv_ratio_50, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(fr_boot_dec.index, inv_ratio_5, inv_ratio_95, color='teal', alpha=0.2)
plt.plot(fr85_boot_dec.index, inv_ratio85_50, color='red', linewidth=4, alpha=0.8, label='RCP 8.5')
plt.fill_between(fr85_boot_dec.index, inv_ratio85_5, inv_ratio85_95, color='red', alpha=0.15)

#label='Theoretical=1',
plt.axhline(y=1, color='blue', linestyle='solid',  linewidth=3, alpha=0.8) #label='Theoretical=1',

#plt.xlabel('Water Year(Apr-Mar)', fontsize=30) #fontweight='bold'
#plt.ylabel('Probabaility Ratio (Log)', fontsize=30)
#plt.title('Chenab drying Annual-Min Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=36) #,fontweight='bold'
plt.yticks(fontsize=36) #,fontweight='bold'
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=18, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=18, width=4)

# Set logarithmic scale for y-axis
plt.yscale('log')
# Set custom y-tick labels for logarithmic scale
yticks_log = [0.2, 0.5, 1, 2, 5, 10]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)

# Increase the thickness of the outer boundary
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Adding a legend
#plt.legend(loc='upper left', prop={'size': 30}) #, 'weight': 'bold'

plt.ylim(0.19, 10)
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:




